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Abstract. 

We assess the suitability of a recent high-resolution central scheme developed bv lKurganov fc Tadmorl l)2000h for 
the solution of the relativistic hydrodynamic equations. The novelty of this approach relies on the absence of 
Riemann solvers in the solution procedure. The computations we present are performed in one and two spatial 
dimensions in Minkowski spacetime. Standard numerical experiments such as shock tubes and the relativistic flat- 
faced step test are performed. As an astrophysical application the article includes two-dimensional simulations 
of the propagation of relativistic jets using both Cartesian and cylindrical coordinates. The simulations reported 
clearly show the capabilities of the numerical scheme to yield satisfactory results, with an accuracy comparable 
to that obtained by the so-called high-resolution shock- capturing schemes based upon Riemann solvers (Godunov- 
type schemes), even well inside the ultrarelativistic regime. Such central scheme can be straightforwardly applied 
to hyperbolic systems of conservation laws for which the characteristic structure is not explicitly known, or in 
cases where the exact solution of the Riemann problem is prohibitively expensive to compute numerically. Finally, 
we present comparisons with results obtained using various Godunov-type schemes as well as with those obtained 
using other high-resolution central schemes which have recently been reported in the literature. 

Key words. Hydrodynamics - Methods: numerical - Relativity - Shock Waves 



1. Introduction 

Relativistic (magneto) hydrodynamical flows are present 
in many astrophysical scenarios involving compact objects 
such as neutron stars or black holes. The production of rel- 
ativistic radio jets in active galactic nuclei, with flow ve- 
locities as large as 99% of the speed of ligh t, involves an ac- 
cretio n disk around a rotating black hole I Begelman et alJ 
Il984l) . The explosive collapse of the core of a massive star 
to a neutron star in type II/Ib/Ic supernovae contains 
fluid moving at re lativistic speeds and strong shocks (see 
e.g. iMullerl l|1998|) and references therein). Scenarios such 
as the collapse of massive stars to black holes in collap- 
sars, or coalescing neutron star binaries have been pro- 
posed as possible candidates to power 7-ray bursts , with 
bulk Lorentz factors larg er than about 100 (see e.g. lPirar] 
l|l999lk lAlov et all l|2000h and references therein). 

A powerful way to improve our understanding of the 
physical mechanisms which operate in those astrophysi- 
cal systems is through (magneto) hydrodynamical rela- 
tivistic simulations. Numerical simulations of relativistic 



flows, both in Minkowski spacetime as in strong gravi- 
tational field scenarios, have received considerable atten- 
tion in recent years (for a comprehensive list of refer- 
ences the reader is a ddre ssed to the recent reviews of 
iMarti & Miiilel lj2003h and lFonn l|200^ h A consensus that 
has slowly emerged is that a class of conservative finite- 
difference schemes based upon Riemann solvers is partic- 
ularly well-suited to solve the hyperbolic system of the 
general relativistic hydrodynamic equations. Such schemes 
are commonly known as Godunov-type schemes, a class 
of the s o-called hig h-resolution shock- capturing (HRSC) 
schemes l|Toroll 19971) . The knowledge of the characteristic 
speeds of the system of equations, i.e., the eigenvalues of 
the Jacobian matrices associated with the fluxes of the 
equations, together with, in most cases, the correspond- 
ing eigenvectors, is the key building block of any Riemann 
solver. A hydrodynamics code must be able, in particu- 
lar, to resolve complex flows in which strong interacting 
shocks could arise. HRSC schemes written in conserva- 
tion form have proven to fulfill the important requirement 
of capturing the profiles of strong discontinuities in a few 
numerical zones without introducing spurious oscillations. 
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The extension of HRSC schemes based on Riemann 
solvers from Newtonian hydrody namics to g enera l 
relativity was first discussed by iMarti et all l|l99l[) 
where one-dimensional test simulations were performed. 
Further contributions towards completely multidimen- 
sional formulations either adopting the 3+1 splitting 
of the spacetime or fully covar i ant formalisms fol- 
lowed llEulderink fc Mellemal Il995l iBanvuls et alJ Il997t 



Font et alhOOOUPapadoDOulos fc Fontl2000l: llbanez et, alJ 
20011) . Nowadays, many of the astrophysical scenar- 
ios mentioned above have been investigated numerically 
using the equations of relativistic hydrodynamics and 
state-of-the-art HRSC sche mes ( s ee, e. g. references in 
IMarti' fc Muliel l|2003j) and iFontl (|2003^ K Furthermore, 
there is an increasing list of relativistic Rieman n solvers 
available in the literature ijMarti fc Mullerl 12003). includ- 
ing t he exact Riemann solver in the specia l relativistic 
limit iMarti fc Mullerlll994 IPons et alJl2000l) . 

The knowledge of the characteristic structure of a hy- 
perbolic system of conservation laws is recommendable, 
irrespective of the particular algorithm to solve it numer- 
ically. From the theoretical point of view such knowledge 
allows to prescribe physically consistent boundary con- 
ditions and helps to understand the physical properties 
of the system. By knowing the characteristic structure of 
the relativistic hydrodynamics equations it was possible 
to find the exact solution to the Riema nn problem, both 
one-di mensional llMarti fc M iillerl fl9 94) and multidimen- 
sional ijPons et al.ll200o(l . As a result, it was possible to 
know how the tangential velocities modify the flow solu- 
tion in relativity in comparison with the Newtonian case 
llPons et a,l.l2000URezzolla fc Zauottil2002HRezzolla et al J 
120031) . 

On the other hand, an alternative approach to up- 
wind methods for solving hyperbolic systems of conserva- 
tion laws by means of non-oscillatory high-order symmet- 
ric total-variation d iminishing (TVD) schemes emerged 
in the mid 1980s 



iNessvahu fc Tadmorl 



Davisl 119841: iRoel Il984t lYeel Il987t 
19901) (see also lYeel lll98fl) and lTorol 



( 1997J) and references therein). Broadly speaking these ap- 
proaches are based either on Lax-Wendroff second-order 
scheme w it h addition al dissipative terms ijDavisI 1984; 
IRoel Il984t lYed Il987|) or on non-oscillatory high-order 
extensions of Lax-Friedric hs first-order central scheme 
l|Nessvahu fc Tadmorl ll990) . One of the nicest properties 
of this approach is that it exploits the conservation form of 
the Lax-Wendroff or Lax-Friedrichs schemes to yield the 
correct propagation speeds of all nonlinear waves appear- 
ing in the solution. Furthermore, this procedure sidesteps 
the use of Riemann solvers, which results in enhanced 
computational efficiency in multidimensional problems. 
Its use is, therefore, not limited to only those systems of 
equations where the characteristic information is explicitly 
known or to systems where the solution of the Riemann 
problem is prohibitively expensive to compute. This ap- 
proach has rapidly developed during the last decade to 
reach a mature status where a number of straightfor- 
ward central schemes of high order can be applied to any 



nonlinear hyperbolic system of conservation laws. In par- 
ticular, the typical results obtained for the Euler equa- 
tions of Newtonian hydrodynamics show a quality com- 
parable to that of HRSC schemes at the expense of a 
small loss o f sharpness of the solution at discontinuities 
l|Tordll99^1 . An up-to-date summary of th e status an d 
applications of this approac h is discussed in iTorol lll997l) : 
iKureanov fc Tadmorl $200$ : iTadmorl l|200l|) . 

In recent years there have been various success- 
ful attempts at applying high-order central schemes to 
solve the relativistic (magneto) hydrodynamics equa- 
tions l | Koide et aT] |l996t iDel Zanna fc Bucciantinil l2002t 
lAnninos fc Fragile! l2003j) . In the context of special and 
general relativis tic magnetohydrodynamics (MHD), K oide 
and coworkers l|Koide et alJll998l iKoide et alJl2002l) ap- 
plied a second-order centra l schem e with nonlinear dissi- 
pation developed by iDavij l|l984|) to the study of black 
hole a ccretion and formation of r elativ istic jets. More re- 
cently IDel Zanna fc Bucciantinil ((20021) assessed a third- 
order convex essentially non-oscillatory central scheme 
in multidimensional special relativistic hydrodynamics, 
later extended to relativistic MHD in l5eT Zanna et alJ 
l|2003|) . These authors obtained results as accurate as 
those of upwind HRSC schemes in standard tests (shock 
tubes, shock reflection test). Yet another central s cheme 
has been lately considered by lAnninos fc Fragile] l(2003l ) 
in one-dimensional special and general relativistic hy- 
drodynamics, where simi l ar res ults to those reported by 
IDel Zanna fc Bucciantinil l|2002l) are discussed. 

The aim of this paper is to assess, as a proof of principle 
and motivated by the recent stir of activity on this topic, 
the validity of a particularly efficient finite-difference cen- 
tral scheme written in conservation form for the solution of 
the relativist ic hydrodynamic equations. T his scheme was 
developed bv lKurganov fc Tadmorl 1 2000l) and has proven 
very accurate for solving different hyperbolic systems of 
conservation laws, including the Newtonian hydrodynam- 
ics equations. To reach our aim we perform a number 
of one- and two-dimensional standard numerical experi- 
ments in flat spacetime, such as shock tubes and t he rela- 
tivist ic version of the so-called flat-faced step test l)Emervl 
1968|). We also present an astrophysical application of our 
adopted central scheme, namely the propagation of a rela- 
tivistic jet. The simulations show the remarkable capabil- 
ities of the numerical scheme to yield satisfactory results, 
comparable to those obtained by HRSC schemes based 
on Riemann solvers, even well inside the ultrarelativistic 
regime. 

The organization of the paper is as follows: in Section[21 
we remind the reader of the form of the system of equa- 
tions of special relativistic hydrodynamics. Section |3 de- 
scribes briefly the numerical scheme we use. Next, in 
Sections01and|niwe present the results of our one- and two- 
dimensional simulations, respectively. SectionEJcontains a 
quantitative comparison between the central scheme and 
various of the most widely used Riemann solver-based 
HRSC schemes available. Finally, the summary of our in- 
vestigation is given in Section In all simulations pre- 
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sented we choose units such that the speed of light c is 
unity. 

2. Relativistic Hydrodynamic Equations 

In Minkowski spacetimc and Cartesian coordinates (t, x l ), 
the local conservation laws describing the motion of a rel- 
ativistic fluid can be cast as a first-order flux-conservative 
system of the form 

3U(w) 3T(w) 

dt dx* [ 1 

(Latin indices run from I to 3.) In the above equations 
U and P are, respectively, the state vector and the flux 
vector along direction x 1 and are defined as 

U(w) = (D,S\t), (2) 
P(w) = {Dv i ,S j v i +p8 ij ,S i -Dv i ). (3) 

Here, 6^ is the Kronecker delta, p is the fluid thermal pres- 
sure related to the rest-mass density p and specific internal 
energy density e via an equation of state, p — p(p, e), and 
v l is the 3-velocity. The definitions of the evolved quanti- 
ties (relativistic densities of mass, momentum and energy) 
in terms of the primitive variables w = (p, Vi , e) are 

D = pW, (4) 
S* = phW 2 v\ (5) 
r = phW 2 -p-D, (6) 

where h is the specific enthalpy, h = I + e + pj p, and W 
is the Lorentz factor satisfying W = u° = 1/ vl — v 2 with 
v 2 = v l Vi. The 3-velocity components are obtained from 
the spatial components of the 4- velocity as v % — u l /u°. 

3. Numerical Scheme 

The time update of system Q from t n to t n+1 is done 
using an algorithm written in conservation form 



At 



72), 



(7) 



where index i labels the numerical cells. The quantities 
Fi-1/2 and Fj + i/2 are the numerical fluxes at the cell in- 
terfaces i— 1/2 and i+1/2, respectively. In Riemann solver- 
based HRSC schemes those numerical fluxes are computed 
by solving a family of local Riemann problems. Central 
sche mes such as Lax-Friedrichs or Richtmyer avoid this. 

In iKureanov fc Tadmoil |20Q0j) , the authors first con- 
struct a fully-discrete central scheme, by building an in- 
termediate mesh of variable cell length, making use of the 
local speed of propagation at each cell interface a i+1 / 2 , de- 
fined by 



max 



(8) 



where p(A) = maxj ( | A, (.A) | ) , Xi(A) being the eigenvalues 
of the Jacobian matrix A = df/dXJ. In addition, super- 
scripts L and R in the above equation stand for the recon- 
structed values of U at the left and right sides of the cor- 
responding numerical cell (i + 1 and i, respectively). These 



are computed from the reconstructed values of the vector 
of primitive variables defined previously. In particular we 
have implemented in our numerical code the third-order 
reconstr uction procedures provide d by both, the PPM 
scheme llColella fc Woodwardll984h and the PHM scheme 
l|Marauinalll999(l . 

In order to avoid the computation of the Jacobian ma- 
trix of system QJ, one can calculate the partial deriva- 
tives of the flux vector with respect to t he state variables 
numer ically (see iLiu Ta.dmorl ^98); l.liang fc Tadmoil 
l|l998^ 1. The construction of the resulting method is not 
simple, and it is only second order accurate in space. Its 
extension to higher spatial order is quite involved and, 
for the equations of relativistic hydrodynamics, many in- 
termediate calculations are necessary due to the root- 
finding procedure needed to recover the primitive vari- 
ables from the conser ved quantities after a time update 
<|Marti fc Mulleill2003h . The final numerical flux depends 
on a i+ i/ 2 and on the partial derivative of the flux vec- 
tor with respec t to Xj, which can be agai n calculated 
numerically (see iKurganov fc Tadmorl ( 2QQ(t) for details) . 
Preliminary results in one-dimensional simple tests in 
Minkowski spacetime are far more inaccurate than the 
ones obtained when using approximate Riemann solvers. 

However, this first version of the scheme admits a semi- 
discrete form, by setting At — > 0. All new variables con- 
structed on the intermediate mesh of the original scheme 
vanish, resulting in a much simpler and robust scheme. 
Now, we can introduce time differencing again by ap- 
plying a standard time discretization such as a Runge- 
Kutta one. In this way, one can obtain a fully-discrete, 
simple, robust, and characteristic-information-free central 
scheme. Hereafter, we will refer to this scheme as Tadmor's 
scheme, with a numerical flux function given by 



i+1/2 



-[f(Uf +1 ) + f(Uf)] 



uf +1 



u 



(9) 

This numerical flux depends only on the local propaga- 
tion speeds aj+1/2 and, due to its simple form, it can be 
implemented and extended to any spatial order straight- 
fo rwardly. The reader is addre ssed to the original article 
of lKurganov fc Tadmorl |2C)00) for a deeper description of 
the numerical scheme. 



4. One-dimensional tests 

4.1. Riemann problems 

A shock tube problem is a particular Riemann problem 
in which the states at both sides of a given interface are 
at rest. The thermodynamical variables of the fluid are 
discontinuous. When the interface is removed the evolu- 
tion results in four constant states separated by three ele- 
mentary waves, a rarefaction, a contact discontinuity and 
a shock wave. The exact solution of the Riemann prob- 
lem in relativistic hydr odynamics for va nishing tangential 
speeds can be found in lMa,rtf fc Mulle! <|l994h . Compared 
to Newtonian hydrodynamics, in relativity the nonlinear 
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Table 1. Values of the PPM reconstruction parameters 
used in some of the simulations presented in Sections ^ 
and |SJ 



Test 


K 


v w 


>7« 


e« 






£ (2) 


Problem 2 


1.0 


5.0 


0.05 


0.1 


0.52 


10.0 


0.5 


Problem 3 


1.0 


50.0 


0.05 


0.1 


0.52 


10.0 


0.5 


Problem 4 


1.0 


50.0 


0.05 


0.1 


0.52 


10.0 


0.5 


Problem 5 


1.0 


5.0 


0.05 


0.1 


0.52 


10.0 


0.0001 


Problem 6 


1.0 


5.0 


0.1 


0.1 


0.52 


10.0 


0.5 


Jet 


1.0 


1.0 


0.05 


0.1 


0.52 


5.0 


0.1 




velocity addition yields a curved velocity profile for the 
rarefaction wave, whereas the Lorentz contraction narrows 
the constant state in between the shock wave and the 
contact discontinuity. These effects become particularly 
strong in the ultrarelativistic regime, especially the lat- 
ter. Shock tube experiments in relativistic hydrodynamics 
have been considered by many authors to calibrate nu- 
merical code s. An up-to-date summ ary of these efforts is 
presented in lMarti fe Mullerl l|2003|) . 

We present simul ations of four R i emann problems: 
cases (a) and (c) in iMarti fc Mullerl (|l994 - with the 
formation of two shocks and two rarefaction waves, re- 
spectively, and the s hock tube problems 1 and 2 of 
IMarti fc Mullerl l)l996|) . for which the solution consists of a 
rarefaction wave moving to the left and a shock wave mov- 
ing to the right, with a contact discontinuity in between. 
The one-dimensional computational domain extends from 
x = to x = 1. At i = the interface is located at x = 0.5. 
We use an ideal fluid equation of state, p = (T-l)pe, with 
r = 4/3 for the first problem and T = 5/3 for the rest. 

We show results obtained using Tadmor's scheme to- 
gether with a third order Runge-Kutta time discretization 
and the spatial cell reconstruction with which we obtain 
the best results, alternativel y, the parabolic reconstructio n 
used in the PPM sch eme l|Colella fc Woodward! l|l984|) : 
IMarti' fc Mullerl lH 99611. or the hyperbolic one used in the 
PHM scheme l|Marqumal l|l994h . The particular PPM pa- 
rameters we choose for the various tests are reported in 
Table Q] (sec the original article by Colella & Woodward 
( 1984) for information on the meaning of these param- 
eters). For an effective comparison with HRSC meth- 
ods based on (approximate) Riemann solvers, we per- 
form the same tests with three of them (HLLE, Roe, 
and Marquina), implem enting them toget her in our code 
in the way explained in lAlov et all lll999l) . The basic al- 
gorithms of eac h of th ese schemes can e.g. be found in 
IMarti' fc Mullerl l|2003f) . We note that while we employ 
throughout the term "Roe" we are actually referring to 
a Roe-type scheme, in the sense that we use arithmetic 
averaging in the numerical flux computation instead of 
Roe-averaging. Finally, in all figures presented in this sec- 
tion the solid lines indicate the exact solution and the 
different symbols (plus signs, crosses, and circles) indicate 
the numerical solution for the pressure, density and veloc- 
ity, respectively. 
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Fig. 1. Tadmor's scheme (upper panel) and HLLE (lower) 
in the relativistic shock tube problem 1 at t = 0.4 using 
CFL=0.5 and PHM spatial reconstruction. Normalized 
profiles of density, pressure and velocity for the computed 
and exact (solid lines) solutions on an equally spaced grid 
of 400 zones. 



4.1.1. Problem 1 

The initial states are pl = 1-0, pl = 1-0, «l = 0.9 (left) 
and pn = 10.0, /5r = 1.0, vr = 0.0 (right). Figure^shows 
the normalized profiles of the pressure, density and veloc- 
ity at time t = 0.4 on an equally spaced grid of 400 zones. 
The top panel corresponds to Tadmor's scheme and the 
bottom panel to the HLLE Riemann solver. In both cases 
we use PHM reconstruction and a CFL number of 0.5. At 
t = 0.4 the left shock wave is located at x = 0.465, the 
contact discontinuity at x = 0.6 and the right shock is 
located at x = 0.76. 

The conservative central scheme we are using captures 
properly the location and propagation speeds of the dif- 
ferent waves, with an accuracy comparable to the HLLE 
Riemann solver (results with other approximate Riemann 
solvers, not shown in the figure, are also similar). In partic- 
ular it is worth mentioning the sharp resolution attained 
at the discontinuities, especially at the shocks which are 
resolved with the same number of points in either scheme, 
thanks to the use of the PHM third-order reconstruction. 
The small amplitude oscillations observed behind the left 
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Fig. 2. Tadmor's scheme (upper panel) and Roe (lower) 
in the relativistic shock tube problem 2 at t — 0.4, using 
CFL=0.5 and PPM spatial reconstruction (see Table ^ 
for the specific values of the PPM parameters used). 
Normalized profiles of density, pressure and velocity for 
the computed and exact (solid lines) solutions on an 
equally spaced grid of 400 zones. 
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Fig. 3. Tadmor's scheme (upper panel) and Marquina's 
flux formula (lower) in the relativistic shock tube problem 
3 at t = 0.35 using CFL=0.5 and PPM spatial recon- 
struction (see Table \I\ior details). Normalized profiles of 
density, pressure and velocity for the computed and ex- 
act (solid lines) solutions on an equally spaced grid of 400 
zones. 



shock entirely disappear when lowering the value of the 
CFL parameter to 0.3. 

4.1.2. Problem 2 

The initial states are pl = 10.0, pl = 1.0, «l = —0.6 (left) 
and pr = 20.0, pr = 10.0, vr = 0.5 (right). Figure shows 
the normalized profiles of the pressure, density and veloc- 
ity at time t = 0.4 on an equally spaced grid of 400 zones. 
The top panel corresponds to Tadmor's scheme and the 
bottom one to Roe's approximate Riemann solver. We use 
PPM cell reconstruction with specific parameters given 
in Table Both schemes capture properly the location 
and propagation speeds of the different waves. It is worth 
pointing out in particular the fine performance of both 
methods at the contact discontinuity, a feature due to the 
use of the PPM reconstruction. 

4.1.3. Problem 3 

The initial states are pl = 13.3, pl = 10.0 (left) and 
Pr = 0.0, pr = 1.0 (right). Figure |3 shows the nor- 
malized profiles of the pressure, density and velocity at 
time t = 0.35 on an equally spaced grid of 400 zones for 



Tadmor's scheme (top panel) and Marquina's flux formula 
(bottom panel). At t = 0.35 the shock wave is located at 
x = 0.79, the contact discontinuity at x = 0.75 and the 
left and right corners of the rarefaction wave are located at 
x = 0.25 and x = 0.56, respectively. The fluid velocity be- 
hind the shock is 0.714. Fig. [3] clearly shows that the loca- 
tion and propagation speeds of the different waves appear- 
ing in the solution are well captured with both schemes. 
Furthermore, the jumps in the different variables, in par- 
ticular the constant state visible in the density in between 
the shock and the contact discontinuity, are also well re- 
solved. 

The grid resolution employed in our simulations al- 
lo ws for a direct c ompa rison with the results reported 
in iMarti fc Mullerl (^396), wh o used a relativistic exten- 
sion of the PPM method l|Colella fc Woodward! Il984[) 
together with the exact relativistic Rieman n solver, as 
well as with those of iDonat et alJ l|l998[L who em- 
ployed a L JlRSC_s^heme based on Marquina's flux for- 
mula l|Donat fc Marauinal H996) and different high-order 
cell-reconstruction schemes. In addition, we can also com- 
pare o ur r esults with those of iDel Zanna fc Bucciantini 
i|2002h and lAnninos fc Fragile! l|2003h obtained with high- 
order central schemes different to the one we use. In all 
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Fig. 4. Tadmor's scheme (upper panel) and HLLE (lower) 
in the relativistic shock tube problem 4 at t = 0.4 using 
CFL=0.4 and PPM spatial reconstruction (see Table ^ 
for details). Normalized profiles of density, pressure and 
velocity for the computed and exact (solid lines) solutions 
on an equally spaced grid of 400 zones. 



cases, the quality of our computed solution is similar to 
that obtained by those authors. In particular it is worth 
mentioning the sharp resolution attained at the disconti- 
nuities, especially at the contact discontinuity, thanks to 
the use of the PPM reconstruction with the specific pa- 
rameters shown in Table ^ As a compari s on th e third- 
order ENO scheme used by iDonat et al.l l|l998l) shows 
more numerical diffusion than Tadmor's scheme when re- 
solving the contact discontinuity and a similar diffusion 
to the one we obtain for the c ase of the shock wave 
(see Fig. 4 in lDonat et al.l (^998)). Similarly, the captur- 
in g of the contact discontinuity wi th the central schemes 
of iDel Zanna fc Bucciantini (2002) (see their Fig. 1) and 
lAnninos fe Fragile! l)2003|) Tsee their Fig. 2) appears also 
more diffused than in our case. 



4.1.4. Problem 4 

The initial states are now j»l = 10 3 ,pl = 1-0 (left) and 
= 10~ 2 ,pr = 1.0 (right). Figure 0] shows the nor- 
malized profiles of the pressure, density and velocity at 
time t = 0.4 on an equidistant grid of 400 zones, obtained 
when using Tadmor's scheme (top) and HLLE (bottom) 
together with the PPM spatial reconstruction (see param- 
eters in Table The flow pattern is similar to that of 



Problem 3 but the relativistic effects make its computa- 
tion much more exigent. The initial pressure jump of five 
orders of magnitude leads to the formation of a thin and 
dense shell bounded by a leading shock front and a trail- 
ing contact discontinuity - a blast wave. The post-shock 
velocity is 0.96 (W ~ 3.5) while the shock speed is 0.986 
(W ~ 6). Resolving the thin density plateau is a demand- 
ing test for any numerical scheme. 

As in Problem 3 we see that the central scheme we 
use gives the correct propagation speeds of the different 
waves, with an accuracy comparable to that of HRSC 
schemes based on Riemann solvers for the same grid res- 
olution. By direct comparis o n of our results with thos e 
reported in iMarti fc Mulled lll996h . IDonat et all <ll998l) 
IDel Zanna fc Buctitmtm'il l|2002() and lAmimosfc Fragile! 
l|2003j) we see that the overall agreement is similar or 
better. In particular, by comparing the maximum value 
obtained for the density in the thin shell, our central 
scheme achieves a 76% of the exact result. The agree- 
ment found between a Godunov-type scheme such as 
HLLE and a central scheme such as Tadmor is remark- 
able. The diffusion observed in the shock wave is essen- 
tially identical to what is found in the previous refer- 
ences, while that of the contact discontinuity is now con- 
siderably smaller due to the steepening step included in 
the PPM routines. With a grid of 400 zones the constant 
state in between the shock and t he constant discontinu - 
ity can not yet be a chieved, as inlMarti fc Mji lerl lll99fi|L 
Donat et all <ll99£ klDel Zanna fc Bucciantini! l|2002h . and 
Anninos fc Fragile! Ij2003|) . In order to get fully converged 



results with Tadmor's scheme one needs to use an equidis- 
tant grid of a bout 1400 zon e s, a sm aller number than what 
is reported in lDonat et alJ l)l998!l . The grid requirements 
can obviously be reduced by employing adaptive mesh re- 
finement in the region around the dense shell. 

4.2. Problem 5: Shock Reflection Test 

In this test an ideal cold fluid (e\ = 0) with velocity V\ hits 
a wall. The fluid is thus compressed and heats up, produc- 
ing a shock which starts to propagate off the wall, leaving 
the fluid behind at rest (v2 = 0). Subscripts 1 and 2 stand 
for the fluid states ahead and behind of the shock, respec- 
tively. The post-shock density is an increasing function of 
the initial flow velocity. The compression ratio a = pij P\ 
satisfies 



1 



i r-i 



£2: 



(10) 



where £2 = W\ — 1. In the Newtonian limit this compres- 
sion ratio is independent of the inflow velocity. On the 
contrary, in the ultrarelativistic limit the density of the 
gas behind the shock is unbounded (a ~ W\). 

In our setup the computational domain covers the in- 
terval [0,1] and the wall is placed at x — 0. We use a com- 
putational grid with 100 zones and an ideal fluid equation 
of state with T = 4/3. As usual in the simulations of this 



Lucas-Serrano et al.: A central scheme for relativistic hydrodynamics 



7 



1 


p/ 1U 

w>n,f ' f * Vv ^?:0OO0O0O0O0O0OOOO0X3OOC 3000000; C"X i 


1 


0.8 


X 




0.6 


- p/10 5 




0.4 


c 


< 


0.2 











illUMIIlJLIIIlHIIHJL UIIIIIJIIIIUIIIII1III1III-IIIII-IIWI-WIWIIIIIIIIILIUIIIIUIIIHJIIIIUIIIIII-I-KIMIIII 











0.2 



0.4 



0.6 



0.8 



1 



Fig. 5. Tadmor's scheme in the relativistic shock reflection 
test using PPM reconstruction (see Table ^ f° r details). 
Normalized profile of the pressure and the density for a 
time when the shock has propagated 0.5 units off the wall. 
The solid line is the exact solution. The symbols indicate 
the numerical solution obtained with a grid of 100 zones 
and CFL=0.4. 



test the specific internal energy of the incoming fluid is set 
to a negligibly small initial value, S\ — 10 -10 . 

Figure |S] shows the normalized profiles of the pres- 
sure and the density at a time when the shock has prop- 
agated 0.5 units off the wall. The profiles shown cor- 
respond to an initial velocity v\ — —0.99999 (Wi ~ 
224). The solid line is the exact solution and the crosses 
and circles indicate the numerical one for the pressure 
and the density, respectively. Tadmor's scheme is capa- 
ble of resolving accurately the shock location and the 
associated huge jump. In comparison with results ob- 
tained with HR S C sch e mes based o n Riemann solvers 
llMarti fc Mullerl Il996t bonat et all ll99St lAlov et al 



1999J) or other central scheme s |DeTzann*a fc Buccianth" 



2002J lAnninos fc Fragilell2003h , die numerical diffusion at 
the shock present in our results is similar. In Fig. 8 of 
iDonat et alJ lll998h the shock is resolved in 5-6 zones, as 



in our central scheme, while in Fig. 2 of iMarti fc Mullerl 
()l996[) it is capt ured in only 3 zones. I n part icular, 
in the results of iDel Zanna fe Bucciantinil l|2002() . who 
only simulate a moderately "cold" gas (p = 0.01) and 
use 250 zones, the corner of the shock wave is not as 
sharply resolved a s in ou r simulation (see also Fig. 9 of 
lAnninos fe Fragile] (j2003)). In addition, the typical over- 
heating error present in the density near the wall is ~ 1%, 
somew hat lower than the one obtained i n IDonat et alJ 
l|l998h and lDel Zanna fc Bucciantinil l|2002l) . We note once 
again that our results are obtained using the PPM recon- 
struction procedure (see Table for details). 

As already show n in previous works 
JPel Zanna fc Bucciantinl 120021: lAnninos fc Fragile! 
120031) . high-order central schemes are then indeed able 
to achieve results comparable to Riemann solver-based 
schemes in the ultrarelativistic limit (we can reach Lorentz 
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Fig. 6. Tadmor's scheme (upper panel) and HLLE (lower) 
in the relativistic propagation of a blast- wave with nonzero 
tangential velocities at t — 0.4 using CFL=0.5 and PPM 
spatial reconstruction (see Table^for details). Normalized 
profiles of density, pressure and velocity for the computed 
and exact (solid lines) solutions on an equally spaced grid 
of 400 zones. 



factors arbitrarily large, e.g. W\ ~ 7070, corresponding to 
v = —0.99999999c). The most important property which 
allows for central schemes such as the one presented 
here to succeed in the ultrarelativistic regime seems to 
be, in retrospect, the conservation form of the scheme, 
something which was not taken into account in the pio- 
neer investigations in (ult ra) relativistic hydrodynamics 
I Norman fc Winkler! Il98^ . 



4.3. Problem 6: Shock Tube Tests with non-zero 
tangential velocities 

All shock tube problems we have considered so far involve 
only the velocity normal to the initial discontinuity. A 
higher level of complexity can be achieved by including 
the two other velocity components v y and v z , which can 
be considered together as v± , the velocity component de- 
finc d on the p e rpend icular plane to the x-axis. As shown 
by IPons et all l(2000[) , where the exact solution of such a 
"multidimensional" Riemann problem was derived, the in- 
troduction of this new variable changes the final result of 
the initial Riemann problem. 

Following the previous reference we have simulated the 
propagation of a relativistic blast wave (see M4. 1.4(1 . The 
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final solution depends on the initial value of the tangen- 
tial velocities, wh ich can be taken as free parameters (see 
IPons et all l|200fjl for details). In our particular case we 
choose Vj_ = 0, = 0.99. We use the same configuration 
as in Section 14.11 for the spatial domain and the equation 
of state. 

In Figure|S]we plot the results obtained using Tadmor's 
scheme and HLLE approximate Riemann solver. As in 
the previous tests we use a third order Runge-Kutta time 
discretization and the PPM scheme for the spatial cell- 
reconstruction (see Table Q] for details on the PPM pa- 
rameters) . The solution plotted in Fig. corresponds to 
time t = 0.40, and it shows the distinctive material shell 
bounded by a contact discontinuity and a shock. The shell 
is now thicker than in the case with no tangential velocities 
(Problem 4) and the density jump is larger. We again find 
that Tadmor's scheme captures successfully the correct 
wave patterns, with a very high accuracy at the discon- 
tinuities, due to the use of the PPM reconstruction. We 
note in particular that the contact discontinuity is simi- 
larly well resolved with both schemes, HLLE and Tadmor. 
It is worth emphasizing that despite the further complex- 
ity introduced by the presence of tangential velocities, the 
correct value for the density of the material shell is nev- 
ertheless computed very accurately. 

5. 2-dimensional Relativistic Hydrodynamics 

The simplicity of Tadmor's scheme allows to implement 
it in multidimensions in a straightforward way, by us ing, 
e.g. the so-called method of lines (see e.g. iTorol l)l997l) ). 
In this framework we have thus considered two differ- 
ent tests and an astrophysical application, namely a two- 
dimensional shock tube, the relativistic version of the so- 
called flat-faced step test, and the propagation of a rel- 
ativistic jet. Our results c an be directly compared with 
those of previous authors llDonat et alJUoiiil iMarauinal 
ll993lPel Zanna fc Buccia,ntiniH2002l). 



5.1. Two-dimensional shock tube test 

In this test the computational domain is a square of side- 
length unity divided in four quadrants of equal size with 
constan t initial states each and outflowing boundary con- 
ditions. lLax fe Liul (^98) studied the time evolution of 
all possible initial configurations within the framework of 
Newtonian hydrodynamics. A relativistic version of their 
particular configuration 12, where the four boundaries de- 
fine two co ntact discontinuities and two sho ck waves, was 
proposed in lDel Zanna fc Bucciantinl |2002). We simulate 
the same configuration in order to compare with their re- 
sults. The initial setup is 



f (p,v x ,v y , P ) NE =(0.1,0.0,0.0,0.01), 
(p,v x ,v y ,p) NW =(0.1,0.99,0.0,1.0), 
(p,v x ,v y ,p) sw =(0.5,0.0,0.0,1.0), 

t (p : v x ,v yiP ) SE =(0.1,0.0,0.99,1.0). 




Fig. 7. Tadmor's scheme with PHM reconstruction in a 

two-dimensional shock tube test. 30 iso-contours of the 

logarithm of the density are shown, t = 0.4, CFL=0.5 and 
l _ j_ _ _j_ 

Ax ~ Ay ~~ 400 ' 



Figure \7\ shows the result we obtain using Tadmor's 
scheme with a 400x400 grid at time t = 0.4. We use 
PHM spatial reconstruction, a third order Runge-Kutta 
time discretization and a CFL number of 0.5. The solution 
shows two curved shock fronts evolving in the upper-right 
(NE) quadrant, and a more complicated wave pattern in 
the low er-left (SW) quadrant. The comp arison with the re- 
sults of iDel Zanna fc Bucciantmll l|2002|) (see their Fig. 6) 
reveals that all expected features in the wave solution are 
also obtained with our central scheme. The capturing of 
the curved shocks in the upper-right quadrant is similarly 
well performed by both schemes. However, the elongated 
structure appearing along the diagonal within the two 
shocks is more pronounced in our res ult (in clos e r agre e- 
ment with the Newtonian results bv lLax fc LiiJ l|l998^ L 
Correspondingly, concerning the wave structure in the 
lower-left quadrant we note that the location and reso- 
lution of the bow shock agrees in both schemes but, how- 
ever, the contact discontinuity and the structure in front 
of the oblique shock differ, being perhaps better resolved 
with our scheme. 



5.2. Relativistic flat- faced step test 

A challenging test for two-dimensional hydrodynamical 
codes is the numerical simulation of a supersonic flow in 
a wind tun nel with a flat- faced step, a test originally in- 
troduced bv lEmervl <[l968) to compare various schemes in 
classical fluid dynamics. The initial configuration of this 
test is as follows: the tunnel is three units long and one 
unit wide. The step is 0.2 units high and it is located 
0.6 units from the left-hand end of the tunnel. A Mach 
3 flow (Newtonian definition) is injected through the left 
end of the tunnel. All the computational domain is ini- 
tially filled with an ideal gas with 7 = 7/5 and constant 
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Fig. 8. Tadmor's scheme with PHM reconstruction in the 
flat-faced step test. 30 iso-contours of the logarithm of 
the density are shown. CFL=0.7 and ^ = = ^j. 
From top to bottom: v x = 0.1,0.9,0.99, and 0.999, and 
t = 55.1,6.0,4.45, and 4.26 

density p{x,y) — 1.4. The only non-vanishing initial ve- 
locity component is the horizontal (x) one, which we leave 
as a free parameter in order to consider different regimes, 
from low Lorentz factors up to ultrarelativistic situations. 
Reflecting boundary conditions are applied along the walls 
of the tunnel as well as on the boundary defined by the 
step. Correspondingly, outflow (zero gradient) boundary 
conditions are used on the right-hand end of the tunnel. 

The corner of the step is a singular point of the flow 
since it is the center of a rarefaction fan. It is well known 
that linearized Riemann solv ers need an entrop y fix in 
the vicinity of the corner (see iDonat et 

ID fjflflffl for de- 

tails) in order to minimize the numerical errors generated 
around it which may affect the entire flow globally. We 
note that central schemes such as the one we use do not 
need any further adjustment to pass this test without di- 
minishing the quality of the results. 

Figure|Hlplots the results obtained by Tadmor's scheme 
in a rectangular grid of 120 x-cells and 40 y-cells, using 
a third order Runge-Kutta time discretization and PHM 
spatial reconstruction. We note that the reconstruction is 
now made on the proper velocity components of the gas 
(which are unbounded), in order to avoid exceeding the 




1 2 



Fig. 9. Tadmor's scheme with PHM reconstruction in the 

fiat-faced step test. 30 iso-contours of the logarithm of 

the density are shown. v x = 0.99, t = 4.5, CFL=0.7 and 
l i _ i 

Ax ~ Ay ~ 80 



speed of light during the recovery of the primitive quanti- 
ties. From top to bottom Fig. [5] shows a snapshot in the 
evolution of the flow for the values of the initial velocity 
v x = 0.1,0.9,0.99, and 0.999. The evolution of the fluid 
and the shock reflections pattern are similar to the one 
found in the Newtonian case, but the bow shock moves 
faster in the relativistic case, eventually leaving the com- 
putational domain, the sooner the larger the inflow veloc- 
ities are. The results obtained with Tadmor's scheme are 
very satisfactory and again comparable to the ones ob- 
tained by Riemann solver-based methods which make use 
of the characteristic information of the syste m of equations 
(see e. g . the correspondi ng figures reported in lDonat et alJ 
ill998}klMarauinal jl999|) L All special features (bow shock, 
rarefactions fan, etc.) are well resolved and correctly cap- 
tured. We note that by refining the grid resolution we 
do not encounter the numerical pathologies commented in 
iMarauinal (^)99) , as we show in Figure |5] where a grid of 
240 x-cells and 80 y-cells was used. It is worth comment- 
ing that when using a seco nd order spa tial reconstruction 
such as the MC limiter fsee lTorol (^97)) and a second or- 
der Runge-Kutta time discretization, the results are still 
acceptable, but with a lower resolution at the discontinu- 
ities. 

5.3. An astrophysical application: propagation of 
relativistic jets 

The simulation of relativistic jets (a field of resea r ch tha t 
started about a decade ago; see e.g. iMarti et alJ 1(19 941 1 
has now reached its maturity, gradually incorporating 
more elaborate ingredients such as three dimensional ef- 
fects, magnetic fields, realistic equations of state, and 
emission processes (see IMarti fc Muiierl d2003h and refer- 
ences therein) . As a sample of an astrophysical application 
of our central scheme, we present in this section numeri- 
cal simulations, using both Cartesian and cylindrical co- 
ordinates, of the propagation of a relativistic jet in two 
spatial dimensions through a homogeneus environment, 
discussing briefly the main results. For the sake of compar- 
ison the initial d ata considered in our simulation s are the 
same as those of iDel Zanna fc Bucciantiml l|2002|) . Hence, 
a light (jet-to-ambient density ratio, 0.01), relativistic (jet 
Lorentz factor, 7.1) jet with internal Mach number 17.9 
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Fig. 10. Tadmor's scheme with PPM reconstruction in a 
relativistic jet simulation. Snapshots of the rest mass den- 
sity distribution (in logarithmic scale) showing, from top 
to bottom, the propagation of the jet at t = 20, 40, 60, 80. 
In the last snapshot, the terminal shock, the contact 
discontinuity and the bow shock are found at z = 
29.0, 31.2, 32.6, respectively. 



(highly supersonic) is injected into a homogeneous, static 
external medium. The jet and the ambient medium are 
in pressure equilibrium. We use an ideal gas equation of 
state with 7 = 5/3 to represent both jet and ambient 
medium. Reflecting boundary conditions are used at the 
jet symmetry axis and a combination of outflow and reflec- 
tion conditions at the boundary above the jet inlet (at the 
left corner of the numerical grid) . Pure outflow boundary 
conditions are imposed in the remaining boundaries. The 
simulations are performed with a resolution of 20 cells per 
jet radius and with a CFL of 0.3. 

Figure EH shows a series of snapshots of the rest mass 
density distribution covering the evolution of the jet up 
to t — 80. This simulation is performed using cylindri- 
cal coordinates (r, z) to discretize the computational do- 
main, which is 45 units long in the z-direction and 25 
units wide in the r-direction. All structural features usu- 
ally appearing in jet simulations are clearly identified in 
this figure. A supersonic beam extends from the jet noz- 
zle to the point of impact with the ambient medium. At 
that point the jet ends up in a complex structure formed 



by a terminal planar shock (Mach shock), a contact dis- 
continuity separating the jet material from the shocked 
ambient medium, and a bow shock. This shock forms due 
to the supersonic propagation of the jet with respect to 
the ambient, allowing for the jet to evolve in a cavity of 
hotter and lighter shocked ambient material. Finally, the 
beam material which stops at the jet end forms a cocoon 
sourrounding the beam. The coccon is separated from the 
shocked ambient medium by a contact discontinuity in 
which Kelvin-Helmholtz instabilities develop. The veloc- 
ity of advance of the jet in the ambient medium is governed 
by the balanc e of momentum at the jet/ambient impact 
region (see e.g. lMarti et al.1 l)l997|) l. For the jet under con- 
sideration a theoretical advance speed of 0.44 is expected 
(to be compared with a mean speed of 0.38 found in the 
simulation) . The lateral expansion of the cavity as well as 
the evolutio n of pressure and density follows the theoreti- 
cal model of lBeeelman fc Ciofil l|l989f) with good accuracy. 

It is worth commenting on the absence of carbun- 
cle ahead of the bow shock in the simulations reported 
in Fig. 1101 The carbuncle is a numerical pathology well 
known from blunt body simulations in supersonic gas 
dynamics, which manifests itself in the form of a small 
unphysical protuberance in front of the bow shock. In 
this regard Tadmor's scheme is able to cleanly resolve 
(carbuncle-free) the leading bow shock in the jet with 
an accuracy comparable to that of more sophisticated 
Godunov-type schemes such as Marquina, and clearly su- 
perior to that of Roe-type approximate solvers, which 
so metimes admit this kind of spurious solution as shown 
in lDonat etafl l|l998h . 

Figure ^2 shows the last computed time of our model 
(t = 100) in 2D cylindrical and planar coordinates. The 
morphological elements found in both cases are similar 
as well as the dynamics governing the jet propagation. 
The change in geometry is the responsible of the different 
aspect of the cavity (more elongated in the case of the 
cylindrical jet). 

6. Quantitative comparison 

In this section we present a quantitative comparison be- 
tween Tadmor's scheme and the HRSC methods based on 
(approximate) Riemann solvers we have used. To this aim 
we first compare the numerical and analytic solutions, and 
show the errors of the various methods for the shock tube 
problem 3 and problem 6 (propagation of a blast wave 
with non-zero tangential velocities). We also show in each 
case the CPU time per numerical cell and iteration (TCI) 
in order to highlight the computational efficiency of our 
central scheme. 

The results of the comparison are shown in Table 
This table reports the L\ and norms for the den- 
sity for Tadmor's scheme, HLLE, and Roe's approximate 
Riemann solvers. Results for both spatial reconstruction 
procedures (PPM and PHM) are also included. The L\ 
norm errors are computed considering the entire compu- 
tational domain, while the errors refer only to the 
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Fig. 11. Tadmor's scheme with PPM reconstruction in a 
relativistic jet simulation. Rest mass density distribution 
(in logarithmic scale) at t — 100 for cylindrical (top panel) 
and planar (bottom panel) jets. 

part of the grid where the solution is smooth (i.e. the rar- 
efaction wave). The largest errors occur in the postshock 
region. From this table we can quantify and confirm the 
quality of the results shown in the corresponding figures, 
i.e. that Tadmor's scheme has an accuracy comparable to 
that of Riemann solvers based HRSC methods. 



Correspondingly, Table displays the TCI for 
Tadmor's scheme and Marquina's flux formula, which al- 
lows us to check their computational efficiency. We note 



Table 2. Errors in the shock tube problems 3 and 6 
for HLLE and Roe's approximate Riemann solvers, and 
Tadmor's central scheme using two different third order 
spatial reconstruction (CFL=0.5 and N x = 400). 



Test 


Scheme 


Li 


(10^) 


Loo 


(io- a ) 






PPM 


PHM 


PPM 


PHM 


Problem 3 


HLLE 


3.37 


2.81 


2.87 


0.67 




Roe 


3.39 


2.95 


3.44 


1.31 




Tadmor 


3.41 


3.41 


2.58 


1.04 


Problem 6 


HLLE 


22.0 


21.7 


3.53 


1.20 




Roe 


21.7 


21.6 


2.63 


1.20 




Tadmor 


25.2 


22.4 


0.92 


2.69 



Table 3. CPU time per numerical cell and iteration for 
Tadmor's central scheme and Marquina's flux formula for 
the different regimes considered in the text, using third 
order schemes for the cell reconstruction and time update. 



TCI (/is) 


Case 


# Zones 


Marquina Tadmor 


ID 


400 x 1 x 1 


17.0 8.0 


1.5D 


400 x 1 x 1 


32.4 8.8 


2D 


120 x 40 x 1 


72.1 22.7 



first that when writing the numerical code we did not 
take special care in optimization issues other than in its 
most quite apparent aspects. Therefore, the numbers dis- 
played in Table |3| have to be taken as approximate num- 
bers for a standard implementation of a hydrodynamics 
code, amenable to be improved under code optimization. 
From this table we can see that Tadmor's scheme, as ex- 
pected, consumes less time in doing the calculations than a 
method based on the characteristic structure of the equa- 
tions, roughly a factor of 2 for one-dimensional problems, 
a factor 4 for problem 6 where tangential flow velocities 
are also involved in the computation, and a factor 3 for 
a two-dimensional problem (namely, the flat-faced step 
test). 



Some additional comments are in order: firstly, the 
TCI for Tadmor's scheme remains essentially unchanged 
when going from ID to 1.5D (from 8.0 /is to 8.8 /xs), 
which indicates the negligible impact of the computation 
of the numerical flux in the overall number of operations in 
the code. Correspondingly, the factor of two difference in 
Marquina's flux formula (from 17.0 /is to 32.4 /is) is con- 
sistent with the large relative weight of the numerical flux 
step in the overall computation in this algorithm (notice, 
in particular, that the numerical viscosity matrix changes 
its dimension from 3 x 3 to 5 x 5). Secondly, when going 
from 1.5D to 2D the TCI in Marquina's scheme increases 
by roughly a factor 2.2, while this factor is of the order of 
2.5 in the case of Tadmor. A factor two can be easily un- 
derstood due to the presence of the additional dimension 
in the code, which implies an additional "sweep" in the 
corresponding spatial dimension, the cell reconstruction 
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Table 4. L\ norm errors of the density and convergence 
rate (r) under grid refinement for problem 4 at t = 0.4. 
Results for three schemes (Marquina, HLLE, and Tadmor) 
and two reconstruction procedures (PPM and PHM) are 
shown. 



Scheme 


N 


Li(10 


2 ) 

J 




r 






PPM 


PHM 


PPM 


PHM 


IVXdlU LIIXJ-CL 


O U 




9^ 9 








±\JyJ 


21.0 


J-O.tJ 


n ^n 

O.OO 


u.ou 




900 


1 fi 1 


1 ^ 1 
±o. ± 


n 


^9 




400 


8.99 


11.4 


0.84 


0.41 




800 


4.22 


7.22 


1.09 


0.66 




1600 


2.22 


3.87 


0.93 


0.90 




3200 


1.04 


2.06 


1.09 


0.91 


HLLE 


^n 


14.8 










i nn 

-LOO 


9n Q 

^O. a 


1 7 n 

-L f .O 


n w 

O. OO 


n 9Q 




9nn 


1 (S 1 


1 ^ 1 
1U. 1 


n 

o.oo 


n n 

o.oo 




/inn 


8 QQ 

o.yy 


11/1 
11.4 


n S/i 


n /i i 




800 


4.22 


7.22 


1.09 


0.66 




1600 


2.22 


3.87 


0.93 


0.90 




3200 


1.04 


2.06 


1.09 


0.91 


Tadmor 


50 


14.7 


23.6 








100 


20.7 


19.1 


-0.49 


0.30 




200 


16.0 


15.2 


0.37 


0.33 




400 


8.94 


11.0 


0.84 


0.47 




800 


4.19 


7.26 


1.09 


0.60 




1600 


2.21 


3.89 


0.92 


0.90 




3200 


1.04 


2.07 


1.09 


0.91 



and the numerical flux computation being also computed 
twice. This factor can be further modified by adding and 
subtracting two less important factors, namely the access 
to memory (which is slower as the arrays are larger) and 
the time update and recovery procedure which are done 
simultaneously for the two spatial dimensions. 

Finally, we show in Table 0] the errors of the density 
under grid refinement using the discrete L\ norm. The 
results reported in this table correspond to problem 4 at 
t = 0.4, and they allow for a comparison of the various 
schemes (Marquina, HLLE, and Tadmor) and reconstruc- 
tion procedures (PPM and PHM). The last two columns 
of the table indicate the convergence rate of the method. 
As the grid is refined it can be seen how the convergence 
rate reaches an order of accuracy of roughly one. This is 
the expected value for problems where discontinuities are 
present. Furthermore, ou r result is also in very go od agree- 
ment with the results of iMarti fc Mullerl (^,996) where a 
relativistic extension of PPM was used in conjunction with 
the exact Riemann solver (see their Table IV). 



Two important conclusions can be drawn from Tabled 
firstly, the transition to the converged solution is faster 
with PPM than with PHM, especially for grid resolu- 
tions finer than 1/400. The corresponding L\ norm er- 
rors are also systematically smaller when using the PPM 
cell-reconstruction routines, roughly a factor of two better 
than PHM when the number of zones is > 400. Secondly, 
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Fig. 12. L\ norm errors of the density under grid refine- 
ment for problem 4 at t = 0.4. Open (respectively, filled) 
symbols correspond to results obtained with PPM (re- 
spectively, PHM). Results for three schemes (Marquina 
(squares), HLLE (circles), and Tadmor (triangles)) are 
shown. For a given grid resolution the accuracy of the 
solution shows no dependence on the numerical scheme 
and a strong dependence on the cell-reconstruction. 

both the convergence rate and the errors are pretty much 
independent of the scheme (Marquina, HLLE, or Tadmor), 
and irrespective of the grid resolution employed. This re- 
sult indicates that for a given scheme (cither upwind or 
central) written in conservation form, it is the reconstruc- 
tion what greatly helps in order to gain accuracy. These 
conclusions can also be inferred from Fig. 1 121 which shows 
the convergence rate in logarithmic scale. It is however 
important to stress that this result applies to the specific 
test we have considered (problem 4), for which the wave 
structure of the solution is particularly simple and the er- 
rors are dominated by the presence of discontinuities. It 
remains to be seen if our conclusion still holds in multi- 
dimensional problems involving more complex flows and 
wave interactions (not necessarily aligned with the com- 
putational grid) as those appearing, for instance, when 
simulating astrophysical jets. 

7. Summary 

In this paper we have assessed the validity of a particular 
finite-differ ence central scheme in cons ervation form de- 
veloped bv iKurganov fc Tadmorl |2000), for the solution 
of the relativistic hydrodynamic equations. The computa- 
tions have been restricted to one and two spatial dimen- 
sions in flat spacetime. Standard numerical experiments 
such as shock tubes, the shock reflection test, and the 
relativistic version of the flat-faced step test have been 
performed. As an astrophysical application of the cen- 
tral scheme we use we have presented two-dimensional 
simulations of the propagation of relativistic jets using 
both Cartesian and cylindrical coordinates. The simula- 
tions have shown the capabilities of Tadmor's scheme to 
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yield satisfactory results, comparable to those obtained by 
HRSC schemes based on Riemann solvers, even well inside 
the ultrarelativistic regime. We have proposed to use high- 
order reconstruction procedures such as those provided 
by the PPM s cheme dColella fc Woodwardlll984|) and the 
PHM scheme (jMarauinalll99^i . This is essential to keep 
the inherent diffusion of central schemes at discontinuities 
at reasonably low levels. 

The novelty of our approa ch (shared by recent ear- 
lier works in the literature l|Del Zanna fc Bucciantini! 
120021 lAnninos fc Fragile! 12003^ relies on the absence of 
Riemann solvers in the solution procedure. Earlier pio- 
neer approaches in the field of relativistic hydrodyn amics 
l|Norman fc Winklerlll986HCentrella fc Wilsor]ll984f) used 
standard finite-difference schemes in conjunction with ar- 
tificial viscosity terms to stabilize the solution across dis- 
continuities. Those approaches, however, only succeeded 
in obtaining accurate results for moderate values of the 
Lorentz factor (W ~ 2). A key feature lacking in those 
earlier investigations was to write the system of equations 
and the numerical scheme in conservation form. To the 
light of our findi ngs and, in addition, to the recent re- 
sults reported by lAnninos fc Fragile] l|2003|) where artifi- 
cial techniques have also been considered in tandem with 
central schemes, it appears that this was the ultimate rea- 
son preventing the extension of the computations to the 
ultrarelativistic regime. 

The use of high-order central schemes for nonlin- 
ear hyperbolic systems of conservation laws is increas- 
ing in recent years since the seminal work in the sec- 
ond half of the 1980s llDavisI Il983 iRoel Il984l lYeel 
Il987t iNessvahu fc Tadmorl Il990f) . Anile and coworkers 
ijAnile et^?Tl200*o|) appliedthe central scheme SLIC in the 
context of the time-dependent hydrodynamical semicon- 
ductor equations obtaining satisfactory results as well for 
a system whose characteristic information is not avail- 
able. In the context of s pecial and general relativistic 
MHD iKoide et all l|l996l ll99S|) applied a second-order 
centra l sche me with nonlinear dissipation developed by 
ID avid l|l984h to the study of relativistic extragalactic jets 
and black hole accretion. This scheme has been lately 
applied by iMizuno et ail l)2003|) in general relativistic 
MHD simulations of the gravitational collapse of a mag- 
netized rotating ma ssive star as a model of gamma ray 
bursts. Also recent ly Del Zanna fc Bucci antini (2002) and 
lAnninos fc Fragile] l|2003j) have assessed two different high- 
order central schemes in relativistic hydrodynamics, ob- 
taining results as accurate as those of upwind HRSC 
schemes in standard tests. 

The theoretical knowledge of the characteristic in- 
formation of any hyperbolic system of equations is the 
key ingredient guiding the construction of HRSC upwind 
schemes. The different Riemann solvers av ailable in the lit- 
eratur e, either exact or approximate (see lMarti fc Mullerl 
(2003) for an up-to-date list in relativistic hydrodynam- 
ics), all use in the solution procedure the eigenvalues (char- 
acteristic speeds) and/or the eigenvectors (characteristic 
fields) of the Jacobian matrices of the system. Such solvers 



provide a minute amount of diffusion across discontinu- 
ities while at the same time capture the jumps in a self- 
consistent way thanks to the implicit use of the Rankine- 
Hugoniot jump conditions (the so-called shock-capturing 
property). Having such information available is also very 
important for the issue of imposing boundary conditions 
in regions where a priori there could exist some ambigu- 
ity. The knowledge of the behavior of the characteristic 
speeds and, therefore, the local directionality of the flow, 
greatly simplifies this task. 

Needless to say, the alternative of using high-order 
central schemes as the one discussed here instead of up- 
wind HRSC schemes becomes apparent when the spec- 
tral decomposition of the hyperbolic system under con- 
sideration is not known. The straightforwardness of a cen- 
tral scheme makes its use very appealing, especially in 
multi-dimensions where computational efficiency is an is- 
sue. Perhaps the most important example in relativistic 
astrophysics is the system of (general) relativistic magne- 
tohydrodynamic equations. Despite some r ecent progress 
has been accomplished in recent ye ars l|Romero et alJ 
Il997t iBalsaralEoOlt lKomissarovlll999|h much more work 
is certainly needed concerning their solution with HRSC 
schemes based upon Riemann solvers. Meanwhile, an ob- 
vious choice is the use of central schemes. 
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